MultiNicheNet Review: Differential connectomics analysis BC dataset
Daniele Pittari
2023-11-25
setwd("./bc_dataset")
library(multinichenetr)
library(Seurat)
library(Connectome)
library(ggplot2)
library(cowplot)
library(ComplexHeatmap)
library(SingleCellExperiment)
library(dplyr)
source('../DifferentialConnectome.wrapper.r', local=T)
dir.create(file.path(getwd(), "result_tables"), showWarnings = FALSE)1 Load the full cohort
In this case, the original raw counts matrix will be utilized (discrete counts).
# Load data sce
sce_subset_bc<-readRDS("D:/Data/datasets/BreastCancer/data/sce_cohort1_updated.rds")
sce_subset_bc<-multinichenetr::alias_to_symbol_SCE(sce_subset_bc, "human") %>%
multinichenetr::makenames_SCE()# Convert to Seurat Obj and assign discrete counts
seurat_subset_sce <- CreateSeuratObject(counts = SingleCellExperiment::counts(sce_subset_bc), meta.data = as.data.frame(colData(sce_subset_bc)))
# free RAM space
rm(sce_subset_bc)
Idents(seurat_subset_sce) <- seurat_subset_sce@meta.data$subType
# Re-assign idents (prevents issues during connectomics analysis)
seurat_subset_sce@meta.data$ident<-Idents(seurat_subset_sce)
table(Idents(seurat_subset_sce))##
## Endothelial_cell Cancer_cell CD8T CD4T
## 6256 53451 18271 19239
## Myeloid_cell NK Fibroblast macrophages
## 6502 2271 22571 6898
## T_cell B_cell CD4REG gdT
## 3093 11500 5798 2152
## Mast_cell DCs monocytes pDC
## 1073 2211 606 20##
## Endothelial_cell Cancer_cell CD8T CD4T
## 6256 53451 18271 19239
## Myeloid_cell NK Fibroblast macrophages
## 6502 2271 22571 6898
## T_cell B_cell CD4REG gdT
## 3093 11500 5798 2152
## Mast_cell DCs monocytes pDC
## 1073 2211 606 20##
## Endothelial_cell Cancer_cell CD8T CD4T
## 2154 20755 3071 3999
## Myeloid_cell NK Fibroblast macrophages
## 2397 488 10754 1764
## T_cell B_cell CD4REG gdT
## 893 2719 1053 251
## Mast_cell DCs monocytes pDC
## 213 559 108 1Last, we remove pDCs: in the downstream analysis of OnNE, the number of pDC is 1, which is not comaptible for allowing the differential connectome analysis.
not_pdc_ids<-row.names(filter(seurat_subset_sce@meta.data, ident!='pDC'))
seurat_subset_sce<-seurat_subset_sce[,not_pdc_ids]
2 Contrasts of interest
The analysis will focused on two aspects that can be highlighted by different contrasts:
PDL1 communication in preE vs PreNE
Chemotactic and (lymph)angiogenic interactions in E (on-pre) vs NE (on-pre)
2.1 PDL1 communication in preE vs PreNE
- CD274 (PDL1) was reported to have higher expression in subpopulations of T cells in expanders patentients
- CD274 interaction was predicted using CellPhoneDB in the contrasts between PreE-PreNE
- crucial interplay IFNG-PDL1 signaling, which rapresents the disease hallmark
2.1.0.1 Differential connectomics analysis principles
Here, we approach the analysis via the application of the standard pipeline provided in the differential connectomics vignette.
In essence, this approach allows to compare two conditions in a single differential connectome (weighted, directed network that describes LR pairs as weighted edges and cell idents as vertices). In this specific case, the authors of Connectome provided the a specific scoring system, called perturbation score. In this case, the prioritization focuses on promoting the ligand receptor pairs for those senders and recievers cell identities that experience the highest absolute value of fold change variation in the test condition vs control, as described below:
In addition, the differential connectomics analysis also computes the the weight norm log fold change metrics. This metrics is the sum of the logarithm of the fold change expression between coondition and control for the ligand and the receptor in the LR pair.
The differential connectomics analysis output dataframe stores these two scores as “score” and “weight.norm.lfc” respectively.
2.1.1 Performing differential connectomics analysis of PreE vs PreNE
# Appliacation of the standard differential connectomics analysis, according to vignette preE-vs preNE
PreE_vs_PreNE<-DifferentialConnectome.wrapper(seurat_subset_sce, "expansion_timepoint", "PreE", "PreNE")
# write results
openxlsx::write.xlsx(arrange(PreE_vs_PreNE, desc(weight.norm.lfc)),
"result_tables/PreE_vs_PreNE_connectome_DB.xlsx", rowNames = FALSE)2.1.1.1 Evaluation of the presence of inferred CD274-mediated in the output
CD274 is missing among the results.
Possible reasons:
Not called DE during DEA for filtering of DifferentialConnectome
Not in Connectome prior knowledge LR DB (FANTOM5)
2.1.3 Differential connectomics analysis using NicheNet v2 DB
# Download LR DB of NicheNetV2
lr_network = readRDS(url("https://zenodo.org/record/7074291/files/lr_network_human_21122021.rds"))
dplyr::filter(lr_network, lr_network$from%in%"CD274" | lr_network$to%in%"CD274")NicheNet V2 DB encodes multiple CD274-related LR pairs.
PreE_vs_PreNE_nichenet<-DifferentialConnectome.wrapper(seurat_subset_sce, "expansion_timepoint", "PreE", "PreNE", db=as.data.frame(lr_network))
# Write xlsx based on desc weight norm lfc, not filtered
openxlsx::write.xlsx(arrange(PreE_vs_PreNE_nichenet, desc(weight.norm.lfc)),
"result_tables/PreE_vs_PreNE_connectome_via_nichenetLR_V2_DB.xlsx", rowNames = FALSE)2.1.4 Evaluation of CD274 pairs and conclusions
# Show CD274 LR reordered by weight.norm.lfc
PreE_vs_PreNE_nichenet%>%
arrange( desc( weight.norm.lfc))%>%
mutate(rank_wnlfc= rank(desc(.$weight.norm.lfc)))%>%
filter(ligand %in%"CD274" | receptor%in%"CD274")%>%
select(6:18)Conclusions: the analysis performed using NicheNet V2 prior model allows to capture and prioritize CD247 LR interactions. However, in the connectome paper, in the context of differential connectomics analysis the authors suggest:
(…)It is generally useful to limit analysis to those edges which have ligand and receptor expression in > 10% of their respective clusters in either the control or the test condition to avoid noisy measurements.; this thresholding functionality is built into CircosDiff and DifferentialScoringPlot.
This thresholding strategy is indeed applied in the mentioned functions as filtering criterias for these downstream visualization steps of the differential connectomics output (not yet applied because beyond of the scope of these markdowns). If we apply this to our final results, all the CD274 related LR interactions get filtered out.
filter(PreE_vs_PreNE_nichenet, pct.source.1 >0.1 & pct.source.2 >0.1 &
pct.target.1 >0.1 & pct.target.2 >0.1 ) %>%
filter(ligand %in%"CD274" | receptor%in%"CD274")2.2 Inference of Chemotactic and (lymph)angiogenic interactions
Connectome does not allow to accommodate an E(On vs pre) vs NE(on vs pre) contrast design directly. Therefore, we will compute and explore the outputs of the E(On vs Pre) and NE(On vs Pre) as workaround. In the following code we perform differential connectomics analysis via:
Differential connectomics analysis of OnE vs PreE and OnNE vs PreNE, standard Connectome vignette
Differential connectomics analysis of OnE vs PreE and OnNE vs PreNE via NicheNet V2 LR DB
# Appliacation of the standard differential connectomics analysis, according to vignette preE-vs preNE
OnE_vs_PreE<-DifferentialConnectome.wrapper(seurat_subset_sce, "expansion_timepoint", "OnE", "PreE")
OnNE_vs_PreNE<-DifferentialConnectome.wrapper(seurat_subset_sce, "expansion_timepoint", "OnNE", "PreNE")
# Appliacation of the standard differential connectomics analysis using NicheNet V2 LR DB
OnE_vs_PreE_nichenet_db<-DifferentialConnectome.wrapper(seurat_subset_sce, "expansion_timepoint", "OnE", "PreE", db=as.data.frame(lr_network))
OnNE_vs_PreNE_nichenet_db<-DifferentialConnectome.wrapper(seurat_subset_sce, "expansion_timepoint", "OnNE", "PreNE", db=as.data.frame(lr_network))
# write results, arranged by weight.norm.lfc
openxlsx::write.xlsx(arrange(OnE_vs_PreE, desc(weight.norm.lfc)),
"result_tables/OnE_vs_PreE_connectome_DB.xlsx", rowNames = FALSE)
# write results, arranged by weight.norm.lfc
openxlsx::write.xlsx(arrange(OnNE_vs_PreNE, desc(weight.norm.lfc)),
"result_tables/OnNE_vs_PreNE_connectome_DB.xlsx", rowNames = FALSE)
openxlsx::write.xlsx(arrange(OnE_vs_PreE_nichenet_db, desc(OnE_vs_PreE_nichenet_db)),
"result_tables/OnE_vs_PreE_connectome_via_nichenetLR_V2_DB.xlsx", rowNames = FALSE)
# write results, arranged by weight.norm.lfc
openxlsx::write.xlsx(arrange(OnNE_vs_PreNE_nichenet_db, desc(OnNE_vs_PreNE_nichenet_db)),
"result_tables/OnNE_vs_PreNE_connectome_via_nichenetLR_V2_DB.xlsx", rowNames = FALSE)2.2.1 Chemotactic interactions
From Bassels et al:
We next identified DEGs between T cells that expand versus T cells that do not expand pre-treatment (Fig. 4a, Extended Data Fig. 5a,b and Supplementary Dataset 5). Non-expanding T cells were more naive (LEF1, SELL, TCF7), while expanding T cells exhibited high effector function (IFNG), immune cell-homing signals (CXCL13, CCL3/4/5), cytotoxicity (GZMB, PRF1, NKG7), antigen presentation (CD74, HLA-DRB1/5, HLA-DQA2)
From MultiNicheNet paper:
Examples of increasing chemotactic interactions are CCL19-CCR7 between fibroblasts and B cells, CXCL12-CXCR4 between endothelial cells and B cells, and CXCL9-CXCR3 between macrophages/fibroblasts/DCs and (regulatory) CD4 T cells
chemotactic<-c('CXCL13', 'CCL3', 'CCL4', 'CCL5', 'CCL19', 'CXCL9' )
# check if genes are encoded in the connectome prior DB
chemotactic%in%connectome.genes## [1] TRUE TRUE TRUE TRUE TRUE TRUEOnE_vs_PreE_nichenet_db%>%arrange( desc( weight.norm.lfc))%>%
mutate(rank_wnlfc= rank(desc(.$weight.norm.lfc)))%>%
filter(ligand%in%chemotactic)2.2.2 Lymphoangiogenic interactions
lymphangiogenesis<-c('VEGFC', 'FLT4', 'ITGA1', 'ITGA2', 'ITGA5')
# check if genes are encoded in the connectome prior DB
lymphangiogenesis%in%connectome.genes## [1] TRUE TRUE TRUE TRUE TRUEOnNE_vs_PreNE%>%
mutate(rank_wnlfc= rank(desc(.$weight.norm.lfc)))%>%
filter(ligand%in%lymphangiogenesis)OnNE_vs_PreNE_nichenet_db%>%arrange( desc( weight.norm.lfc))%>%
mutate(rank_wnlfc= rank(desc(.$weight.norm.lfc)))%>%
filter(ligand%in%lymphangiogenesis)Conclusions: no lymphoangiogenic interactions using NicheNet and mentioned in the paper were found using Connectome.
## R version 4.3.1 (2023-06-16 ucrt)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 19044)
##
## Matrix products: default
##
##
## locale:
## [1] LC_COLLATE=English_Belgium.utf8 LC_CTYPE=English_Belgium.utf8
## [3] LC_MONETARY=English_Belgium.utf8 LC_NUMERIC=C
## [5] LC_TIME=English_Belgium.utf8
##
## time zone: Europe/Brussels
## tzcode source: internal
##
## attached base packages:
## [1] stats4 grid stats graphics grDevices utils datasets
## [8] methods base
##
## other attached packages:
## [1] gtools_3.9.4 plotrix_3.8-2
## [3] tibble_3.2.1 dplyr_1.1.3
## [5] SingleCellExperiment_1.22.0 SummarizedExperiment_1.30.2
## [7] Biobase_2.60.0 GenomicRanges_1.52.0
## [9] GenomeInfoDb_1.36.4 IRanges_2.34.1
## [11] S4Vectors_0.38.2 BiocGenerics_0.46.0
## [13] MatrixGenerics_1.12.3 matrixStats_1.0.0
## [15] ComplexHeatmap_2.16.0 cowplot_1.1.1
## [17] ggplot2_3.4.3 Connectome_1.0.1
## [19] SeuratObject_4.1.4 Seurat_4.4.0
## [21] multinichenetr_1.0.3
##
## loaded via a namespace (and not attached):
## [1] progress_1.2.2 nnet_7.3-19
## [3] locfdr_1.1-8 goftest_1.2-3
## [5] Biostrings_2.68.1 vctrs_0.6.3
## [7] spatstat.random_3.1-6 digest_0.6.33
## [9] png_0.1-8 shape_1.4.6
## [11] proxy_0.4-27 ggrepel_0.9.3
## [13] deldir_1.0-9 parallelly_1.36.0
## [15] MASS_7.3-60 reshape2_1.4.4
## [17] httpuv_1.6.11 foreach_1.5.2
## [19] withr_2.5.1 xfun_0.40
## [21] ggpubr_0.6.0 ellipsis_0.3.2
## [23] survival_3.5-7 memoise_2.0.1
## [25] ggbeeswarm_0.7.2 emmeans_1.8.8
## [27] zoo_1.8-12 GlobalOptions_0.1.2
## [29] pbapply_1.7-2 Formula_1.2-5
## [31] prettyunits_1.2.0 KEGGREST_1.40.1
## [33] promises_1.2.1 httr_1.4.7
## [35] rstatix_0.7.2 globals_0.16.2
## [37] fitdistrplus_1.1-11 rstudioapi_0.15.0
## [39] miniUI_0.1.1.1 generics_0.1.3
## [41] base64enc_0.1-3 zlibbioc_1.46.0
## [43] ScaledMatrix_1.8.1 ggraph_2.1.0
## [45] polyclip_1.10-6 randomForest_4.7-1.1
## [47] GenomeInfoDbData_1.2.10 xtable_1.8-4
## [49] stringr_1.5.0 doParallel_1.0.17
## [51] evaluate_0.22 S4Arrays_1.0.6
## [53] hms_1.1.3 irlba_2.3.5.1
## [55] colorspace_2.1-0 visNetwork_2.1.2
## [57] ROCR_1.0-11 reticulate_1.32.0.9002
## [59] spatstat.data_3.0-1 magrittr_2.0.3
## [61] lmtest_0.9-40 readr_2.1.4
## [63] later_1.3.1 viridis_0.6.4
## [65] lattice_0.21-9 spatstat.geom_3.2-5
## [67] future.apply_1.11.0 genefilter_1.82.1
## [69] XML_3.99-0.14 scattermore_1.2
## [71] scuttle_1.10.2 shadowtext_0.1.2
## [73] RcppAnnoy_0.0.21 class_7.3-22
## [75] Hmisc_5.1-1 pillar_1.9.0
## [77] nlme_3.1-163 iterators_1.0.14
## [79] caTools_1.18.2 compiler_4.3.1
## [81] beachmat_2.16.0 stringi_1.7.12
## [83] gower_1.0.1 tensor_1.5
## [85] minqa_1.2.6 lubridate_1.9.3
## [87] plyr_1.8.9 crayon_1.5.2
## [89] abind_1.4-5 scater_1.28.0
## [91] blme_1.0-5 locfit_1.5-9.8
## [93] sp_2.1-0 bit_4.0.5
## [95] graphlayouts_1.0.1 UpSetR_1.4.0
## [97] codetools_0.2-19 recipes_1.0.8
## [99] BiocSingular_1.16.0 bslib_0.5.1
## [101] e1071_1.7-13 GetoptLong_1.0.5
## [103] plotly_4.10.2 remaCor_0.0.16
## [105] mime_0.12 splines_4.3.1
## [107] circlize_0.4.16 Rcpp_1.0.11
## [109] sparseMatrixStats_1.12.2 blob_1.2.4
## [111] knitr_1.44 utf8_1.2.3
## [113] clue_0.3-65 lme4_1.1-34
## [115] listenv_0.9.0 checkmate_2.2.0
## [117] DelayedMatrixStats_1.22.6 Rdpack_2.5
## [119] openxlsx_4.2.5.2 ggsignif_0.6.4
## [121] estimability_1.4.1 Matrix_1.6-1.1
## [123] statmod_1.5.0 tzdb_0.4.0
## [125] tweenr_2.0.2 pkgconfig_2.0.3
## [127] tools_4.3.1 cachem_1.0.8
## [129] RSQLite_2.3.1 RhpcBLASctl_0.23-42
## [131] rbibutils_2.2.15 DBI_1.1.3
## [133] viridisLite_0.4.2 numDeriv_2016.8-1.1
## [135] fastmap_1.1.1 rmarkdown_2.25
## [137] scales_1.2.1 ica_1.0-3
## [139] nichenetr_2.0.4 broom_1.0.5
## [141] sass_0.4.7 patchwork_1.1.3
## [143] coda_0.19-4 carData_3.0-5
## [145] RANN_2.6.1 rpart_4.1.21
## [147] farver_2.1.1 aod_1.3.2
## [149] tidygraph_1.2.3 mgcv_1.9-0
## [151] yaml_2.3.7 DiagrammeR_1.0.10
## [153] foreign_0.8-85 cli_3.6.1
## [155] purrr_1.0.2 leiden_0.4.3
## [157] lifecycle_1.0.3 caret_6.0-94
## [159] uwot_0.1.16 glmmTMB_1.1.8
## [161] mvtnorm_1.2-3 bluster_1.10.0
## [163] lava_1.7.2.1 backports_1.4.1
## [165] annotate_1.78.0 BiocParallel_1.34.2
## [167] timechange_0.2.0 gtable_0.3.4
## [169] rjson_0.2.21 ggridges_0.5.4
## [171] progressr_0.14.0 parallel_4.3.1
## [173] pROC_1.18.4 limma_3.56.2
## [175] jsonlite_1.8.7 edgeR_3.42.4
## [177] bitops_1.0-7 bit64_4.0.5
## [179] Rtsne_0.16 spatstat.utils_3.0-3
## [181] BiocNeighbors_1.18.0 zip_2.3.0
## [183] muscat_1.14.0 jquerylib_0.1.4
## [185] metapod_1.8.0 dqrng_0.3.1
## [187] pbkrtest_0.5.2 timeDate_4022.108
## [189] lazyeval_0.2.2 shiny_1.7.5
## [191] htmltools_0.5.6 sctransform_0.4.0
## [193] glue_1.6.2 factoextra_1.0.7
## [195] XVector_0.40.0 RCurl_1.98-1.12
## [197] scran_1.28.2 gridExtra_2.3
## [199] EnvStats_2.8.1 boot_1.3-28.1
## [201] igraph_1.5.1 variancePartition_1.30.2
## [203] TMB_1.9.6 R6_2.5.1
## [205] sva_3.48.0 tidyr_1.3.0
## [207] DESeq2_1.40.2 gplots_3.1.3
## [209] fdrtool_1.2.17 cluster_2.1.4
## [211] ipred_0.9-14 nloptr_2.0.3
## [213] DelayedArray_0.26.7 tidyselect_1.2.0
## [215] vipor_0.4.5 htmlTable_2.4.1
## [217] ggforce_0.4.1 car_3.1-2
## [219] AnnotationDbi_1.62.2 future_1.33.0
## [221] ModelMetrics_1.2.2.2 rsvd_1.0.5
## [223] munsell_0.5.0 KernSmooth_2.23-22
## [225] data.table_1.14.8 htmlwidgets_1.6.2
## [227] RColorBrewer_1.1-3 rlang_1.1.1
## [229] spatstat.sparse_3.0-2 spatstat.explore_3.2-3
## [231] lmerTest_3.1-3 ggnewscale_0.4.9
## [233] fansi_1.0.4 hardhat_1.3.0
## [235] beeswarm_0.4.0 prodlim_2023.08.28